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Title  Slide:  Introduction  to  Wavelets 

This  half-day  course  on  wavelets  is  intended  to  be  a  general  introduction  to  this  exciting 
new  topic  in  signal  processing.  Wavelets  have  been  proven  effective  for  a  number  of  signal 
processing  tasks,  including  compression  and  edge  and  transient  detection.  The  computational 
complexity  of  a  wavelet  decomposition  is  0(N),  versus  0(/V  log  N)  for  the  FFT,  making  it 
very  competitive.  In  addition,  it  has  several  specific  advantages  over  the  FFT  in  many 
applications.  Two  of  these  applications  will  be  discussed  in  this  course:  transient  detection 
and  image  compression. 

Wavelets  have  been  applied  for  a  wide  variety  of  tasks,  and  MITRE  has  been  involved  in 
wavelet  research  for  a  number  of  years,  leading  to  over  thirty  publications.  This  course  was 
developed  in  conjunction  with  W090  and  W096  staff  (who  ran  a  successful  3-day  wavelet 
course  in  Washington),  and  taught  in  Bedford  by  P.  Topiwala.  This  introductory  course  will 
be  followed  by  a  second  course  covering  the  wavelet  research  and  applications  currently 
being  developed  at  MITRE. 


Course  Outline  1 


•  Preliminaries 

•  Theory 

•  Applications 
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Course  Outline 

This  course  will  be  divided  into  an  introduction  to  the  theory  and  to  the  applications. 
Because  of  the  wide  audience  of  this  course,  a  general  summary  of  preliminary  material  is 
provided  as  well. 
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Preliminaries 


•  Linear  Algebra 

•  Linear  Analysis 

•  Fourier  Transform 

MITRE 

Preliminaries 

The  preliminaries  include  linear  algebra  in  finite  dimensions,  linear  analysis  in  infinite 
dimensions,  and  some  well-known  facts  about  the  Fourier  transform. 
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Linear  Algebra  1 


•  Vector  space  Vover  R,C  v.weV,  a,  pe  R,C=»  av  +  pweV 

•  Finite  dimensions  :  Rn,Cn 

n 

•  Inner  product  :  (v,w)=  X  viwi 

i 

•  Norm:  |v|= 
n  " 

•  Basis  {eih  :  v  =  Xaiei  ,oi«  R,C 

i 

•  Orthonormal  (CM)  Basis:  ai  =(v,»i) 

That  Is,  v=  x(v,ei)ei 

i 

•  Example:  «k  =  (0,»  ,1,-0)  1  inthekth  spot  (standard  basis) 

MITRE 

Linear  Algebra  1 

This  slide  covers  the  basic  notions  of  a  vector  space  and  an  orthonormal  basis,  including 
the  definition  of  an  inner  product,  and  the  standard  orthonormal  basis. 


Linear  Algebra  2 


•  Linear  Map  :  A: V  -*  W 

A(a  Vi  +  p  V2)  =  o  A(vi)  +  P  A(v2> 

•  DualBaais:  V,{ek}  eJ:V-»R,C 


ak(ai)s  8u  *  |q 


•  Matrix  Rep  .  (V,{ek}),  (W,{||}) 
A(ak)®  laufi 


•  Example  (standard)  ek=  el  s(0,-">V--0),  1 1n  the  kth  spot 


A«(ai|)  If  A(ei)=  £■«•! 


/MITRE 


Linear  Algebra  2 

A  linear  map  between  two  vector  spaces  can  be  represented  as  a  matrix,  given  bases  for 
the  vector  spaces. 
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Linear  Algebra  3 


•  V©  W  has  basis  ei,— ,#n,fV--.fm;  dim  =n  +  m 

•  V®Whas  basis  {ti®f|},  1=  1,—  ,n; J=  1,— ,m.  dim«nm 
space  of  matrices  of  size  nxm ,  basis  {eT  ©  f  |} 

•  If  A:V  -»  V,  B:W-»W,  then 

A©B:V©W->V©W  ,  A  ®  B:V  ®  W  -»  V  ®  W 

A  ©B  (v,w)  =  (Av.Bv)  ,  A  ®  B  (v ®  w)  =  Av  ®  Bw 

,  more  involved 

MITRE 

Linear  Algebra  3 

Given  two  vector  spaces  with  bases,  new  vector  spaces  can  be  constructed  using  the 
direct  sum  (0)  or  tensor  product  (®). 
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Linear  Algebra  4 


•  Isomorphism  :  A:  V  -*  W  such  that  thers  exists 
B:W  -»  V  with  BA  «1»,AB=1. 

•  For  matrix  ,  nssd  dat  A  #0.  Matrixes  must  be  square 
GLn(R) ,  GLn(C) 

•  Isometry  or  unitary  map  :  (Avi,Av2  )  *  (vi,V2), 
for  any  v  i,V2«V 

•  For  matrix  ,  need  A*A<  AA+  *  I, 
where  A  *  =  (AT) ,  Hermitlan  conjugate 
Ofn).  U(n) 


V 


MITRE 


Linear  Algebra  4 

Two  vector  spaces  are  equivalent  (or  isomorphic)  if  maps  can  be  found  in  either 
direction,  such  that  their  composition  is  the  identity  map  both  ways.  An  important  related 
concept  is  that  of  a  unitary  map,  which  preserves  inner  products.  This  is  automatically  an 
isomorphism. 
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Linear  Algebra  5 


Hermltian  :  A*  A*  aij  ci{i 

Then  Bt  =  Is  unitary  ,  fort  eR,  and  lirn  Bt  =1 

t-»o 

Diagonalization  :  A  Hermltian  map  A  :V  -+  V  has  an 
ON  basis  {#k}  »  A(ek)  =  «*kek. 

r«i  0  ^ 


A  «-* 
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Linear  Algebra  5 

A  Hermitian  matrix  can  be  diagonalized.  Its  eigenvectors  can  be  chosen  to  be 
orthonormal,  and  its  eigenvalues  are  real. 
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Linear  Algebra  6 

*  Frames:  {©$ 

There  are  numbers  0  <  A  <  B  <  <» 
such  that  for  any  v  in  V, 

A||v||2<I|<v,e^|2<B||v||2 

*  Tight  if  A  s  B;  Exact  if  we  cannot  remove 
any  vector  without  losing  completeness 

(1,0),  (1,1)  an  exact  frame 

(1,0),  (1,0),  (0,1),  (0,1)  is  a  tight  frame,  A=B=2 

*  Tight  exact  frame  with  A=B=1  is  an  ON  Basis 

MITRE 
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Linear  Algebra  6 

While  we  will  mainly  deal  with  orthonormal  bases,  there  is  a  generalization  of  the  notion 
of  bases  that  is  useful:  frames.  A  frame  is  a  set  of  vectors,  often  too  many  to  be  a  basis,  such 
that  the  size  of  any  vector  in  the  space  can  be  bounded  above  and  below  in  turns  of  the  inner 
products  of  that  vector  with  the  vectors  in  the  frame.  A  tight  frame  is  a  basis,  while  an  exact 
frame  is  one  where  the  bound  constants  are  equal.  An  orthonormal  basis  is  a  tight  exact 
frame  with  bound  constants  equal  to  1. 
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Linear  Analysis  1 


•  Infinite  dimensions 


:  I2  =  |(si)|liei|  <-|  complex 


sequences 


•  L3(R)  =  tf:R  -»  c|  J  HI*  <  —K  complex  funetlone 

•  Also LJflO,1D.L2(Rn),La(X>,  XcR" 


•  Example  ,  X  *  S1.  circle  In  R  2 


•  All  these  spaces  ere  isomorphic,  In  fact  unltarlly! 

That  la  because  they  all  have  an  orthonormal  basis. 
Different  "realizations'’  of  one  space  H,  Hilbert  space. 


MITRE 


Linear  Analysis  1 

A  finite  dimensional  space  is  equal  to  the  space  of  n  arbitrary  numbers  (n  =  dimension). 

In  infinite  dimensions,  instead  of  dealing  with  spaces  of  arbitrary  numbers,  the  notion  of 
convergence  enters.  One  has  the  space  of  infinite  sequences  where  the  sum  of  the  norm- 
squared  terms  is  finite.  While  there  are  several  other  natural  generalizations  (e.g.,  the 
functions  on  the  real  line  whose  norm-squared  integral  is  finite),  all  of  these  spaces  are 
isomorphic  to  each  other  (by  a  unitary  map).  That  is  because  they  all  have  orthonormal  bases 
of  the  same  cardinality. 
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r  \ 


Linear  Analysis  2 


• 

La([0,2])  =  La([0,1]u  [1,2])=  La([0,1])©  La([1,2]) 

=  H©  H  =  H! 

• 

La(R)=Laf  U  [n,n+1]l=  ©  La([n,n+1D=  ©Hn  sH! 

^n«-«  J  nm-m  — 

• 

La(R")=La(R)®  -  ®La(R)  =  Hi®---® Hn  =H! 

Although  all  of  these  spaces  are  really  the  same,  their 
dlfferent  aspects  give  rise  to  different  bases. 

MTTRE 

Linear  Analysis  2 

If  we  denote  the  common  space  by  H,  the  fact  that  the  various  definitions  lead  to  the  same 
space  implies  some  interesting  things  about  H  under  direct  sum  (®)  and  tensor  product  (0). 
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Linear  Analysis  3 


•  I,  has  standard  basis  {ek},ak«fo,—  ,1,— ,0,— j 

•  La(R)  has  no  natural  ON  basis! 

Traditionally,  Important  ON  bases  arose  In  relation 
to  natural  Hermltian  operators . 


•  "Particle  on  a  loop  "  L2(s1) 

V2f«Xf,  Eigenvectors  of  Laplace  operator 
{o2**™1 },  n  e  Z  (Integers) 

•  l2(R) — i->l2  ON  basis  {en  }  of  L2  (R) 

f  *-»  {<f»en )} 


MURE 


Linear  Analysis  3 

The  sequence  space  has  a  standard  basis,  but  the  other  “realizations”  of  H  have  no  natural 
bases.  Orthonormal  bases  frequently  arise  in  relation  to  Hermitian  operators.  The  special 
functions  of  mathematical  physics  are  instances  of  this.  In  particular,  the  Fourier  series  basis 
derives  from  the  Laplace  operator  on  a  circle  (or  an  interval).  Once  an  orthonormal  basis  is 
secured,  the  equivalence  to  the  sequence  space  is  constructed  by  taking  dot  products  with 
basis  vectors. 
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Linear  Analysis  4 


•  ViewHasL*(R) 

•  Linear  map  T :L2  (R)  -> L2(R)  can  be  represented  as  an  integral  operator 

(Tf)(  x)  =  J  k(  x,  y)f  ( y)dy, 

some  "function"  k(x,y),  called  the  kernel,  on  R2 

•  The  "function”  k(x,y)  may  have  6*  function  singularities 

•  Example:  T  =  1.  Then  k(x,y)=  S(x-y) 

•  "Resolution  of  Id" :  If  {e„ (x)}  Is  an  ON  basis, 
then  for  identify  map  I ,  k(x  -  y)  -  £en(x)e„(y) 

n 

[Physics:  1s£|n><n|l 

.  ^  "  *  MITRE _ 

Linear  Analysis  4 

In  fact,  every  linear  map  on  L2(R)  can  be  represented  as  an  integral  operator  (although 
the  “kernel”  may  have  delta  function  singularities,  e.g.,  the  case  of  the  identity  operator). 
This  is  the  function  space  version  of  the  fact  that  any  linear  map  can  be  represented  by  a 
matrix. 
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Linear  Analysis  5 

• 

Hermltian:  k(x,y)  =  k(y,x) 

• 

LTI:  k(x,y)  =  k(x  -  y),  some  k:R  -> C 

• 

Convolution:  (f*g)(x)«  J  f(x-  y)  g(y)dy 

So  LTI  systems  are  convolution  operators, 
h  =  Impulse  response  function 

• 

Matrlx:ON  basis  {ek},k  =  l,2," 

T:H-»H  has(ai]),Tei  =  £ai|ej 

1 

MITRE 

Linear  Analysis  5 

In  the  case  of  a  linear  “time-invariant”  operator,  the  kernel  is  really  a  function  of  one 
variable.  Such  an  operator  is  thus  a  convolution  with  a  fixed  function,  called  the  impulse 
response  function  (the  output  of  a  delta  function  input). 


14 


Fourier  Transform  1 


•  Unitary  Transform :  T:H-*H 

(Tv,Tw)=  <v,w>,  ail  v,w  e  H 

Examples:  (T«f)(x)  =  f(x-a)  translation 
(E«f)(x)  =  e2*1"  f(x)  modulation 
(D*f)(x)  =  VS  f(ax>  dilation 

A  nxn  matrix,  det  A  s  1  (SLp  (R» 

Ta:L2  (Rb)-»  L2  (Rn) 

(T*f)(x)=f(Ax) 

MITRE 

Fourier  Transform  1 

Important  examples  of  unitary  transforms  of  the  function  space  L2(R)  are  translation, 
modulation  and  dilation. 
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Fourier  Transform  2 


\ 


•  Fourier  Transform: 

(F  /)(«>)  =  /(<»)=  J  r2**®/  (X )  &r  =  u,  Efl)  d)> 

•  "Miracle"  D/(a>)  =  (2nioi)f<a) 

Turns  differential  equations  in  algebraic  equations 

•  (£^j~  to ^  Parseval's  Thm 

•  /£  Convolution  Thm 


V 


MITRE 
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Fourier  Transform  2 

From  the  definition  of  the  Fourier  transform,  it  is  easy  to  see  that  it  is  unitary  (Parseval’s 
theorem).  An  important  fact  about  the  Fourier  transform  is  that  it  turns  differentiation 
operators  into  multiplication  operators,  simplifying  the  solution  of  differential  equations.  For 
signal  analysis,  the  crucial  fact  is  that  the  Fourier  transform  turns  convolution  into 
multiplication.  These  important  facts  justify  the  central  role  of  the  Fourier  transform  in 
analysis.  We  will  see  that  certain  practical  considerations  will  lead  us  elsewhere,  however. 
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Theory 

•  Gabor  Wavelets 

•  Affine  Wavelets 

•  Time-Frequency  Localization 

MITRE 

Theory 

In  this  section  we  outline  the  theory  of  the  wavelet  transform  as  it  exists  today.  There  are 
basically  two  types  of  wavelets:  Gabor  wavelets,  which  relate  to  translations  and 
modulations;  and  affine  wavelets,  which  deal  with  translation  and  dilations.  There  is  an 
essential  difference  between  these  two  types  of  wavelets  in  the  way  they  divide  the  time- 
frequency  plane.  This  leads  us  into  some  general  considerations  of  localization  in  the  time- 
frequency  plane. 
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Overview 

Wavelets  relate  to  a  great  number  of  topics  in  signal  processing.  In  this  course,  we  will 
cover  the  time-frequency  aspects  of  wavelets,  along  with  considerations  of  sampling, 
filtering,  and  bases  and  frames. 
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History 


\ 


1910  A.  Haar 

1920  J.L.  Walsh 

1946  D.  Gabor 

1977  D.  Esteban,  C.  Galand 

1980-5  P.  Burt,  E.  Adelson 

M.J.T.  Smith,  T.P.  Barnwell 
M.  Vetted! 

1985  A.  Grossmann,  J.  Moriet 
1987  I.  Daubechies,  Y.  Meyer,  S.  Mallat 


V 
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History 

In  1910,  Haar  discovered  an  orthonormal  basis  constructed  from  a  single  function,  subject 
to  discrete  translations  and  dilations.  Another  such  orthonormal  basis  was  not  discovered 
until  1985,  by  Daubechies  and  others.  Meanwhile,  in  1946,  D.  Gabor  laid  the  foundations  of 
the  time-frequency  approach  to  bases,  leading  to  the  short-time  Fourier  transform.  In  1977, 
Estaban  and  Galand,  researchers  at  IBM-France,  devised  a  filtering  method  into  low  and  high 
passbands  that  were  free  of  aliasing  on  reconstruction.  Although  an  FIR  filter  could  not  be 
found  to  satisfy  their  criteria  (a  problem  solved  by  Smith  and  Barnwell  in  1980),  they  laid  the 
foundations  of  the  present  theory  of  wavelets.  Grossmann  and  Moriet  constructed  affine 
wavelet  bases  for  good  signal  decomposition  in  the  context  of  oil  exploration.  Finally,  the 
mathematicians  Daubechies,  Meyer  and  Mallat  constructed  an  elegant  theory  for  this  time- 
scale  approach  to  bases. 
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The  Fourier  Transform 


1 

ft 

* 

4 

1 

1 

1 

L. _ 

1 - 

4 _ : 

Mapoitud*  ot  Fourif  Tiflwin 


•  Fourier  theory  is  well-established 

—  Decomposition  of  signal  into  frequency  components 
(Fourier  coefficients) 

-  Basic  signal  components  given  by  functions  e2,r‘f  * 

MITRE  _ _ 


The  Fourier  Transform 

To  introduce  the  wavelet  basis,  recall  that  the  Fourier  transform  decomposes  a  signal  into 
its  harmonic  components.  This  is  useful  when  the  signal  has  well-defined  periodicities,  such 
as  the  sum  of  two  sine  waves  as  depicted. 
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The  Fourier  Transform  (Concluded) 

In  the  context  of  stochastic  processes,  one  requires  the  assumption  of  stationarity.  While 
the  signal  may  be  random,  statistically  it  has  the  same  harmonies  at  all  times.  For  signals  that 
do  not  satisfy  this  assumption  (no  actual  signal  is  stationary),  so  that  the  harmonic  content  is 
time-varying,  the  Fourier  method  is  inadequate  in  that  it  does  not  register  the  time  structure 
of  harmonics.  As  an  example,  the  Fourier  transform  of  a  musical  piece  would  tell  you  only 
what  notes  occurred,  not  when  they  occurred. 

To  get  around  this  problem,  several  methods  have  been  developed.  Gabor  thought  to 
apply  a  running  “window”  to  the  signal  before  doing  the  Fourier  transform,  now  called  the 
short-time  Fourier  transform.  More  generally,  there  are  ways  to  represent  a  running  spectral 
content  of  a  signal,  not  dependent  on  a  windowing  function.  These  fall  under  the  rubric  of 
time-frequency  distributions,  which  we  will  discuss  later.  Finally,  there  is  the  time-scale 
approach  of  the  affine  wavelets,  which  is  the  focus  of  this  course. 


21 


Non-Stationary  Signal  Processing 


Time-dependent  spectrum 


SpMch  Signal 

•  Wavelets  are  a  "new”  tool  for  time-frequency  analysis 

•  Two  classes  of  wavelets 

-  Gabor  wavelets:  windowed  Fourier  transform 

—  Affine  wavelets:  provide  multiscale  decomposition 
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Non-Stationary 

These  new  wavelet  tools  allow  the  analysis  of  nonstationary  signals,  which  abound  in 
nature.  A  classic  example  is  speech  (or  music),  whose  “information”  is  clearly  in  the 
changing  spectral  content 
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Gabor  Wavelets 


/ 


N 


•  r  1  i  4  i 


Ttme-Frequtncy  Plan* 


•  Window  functions  have  fixed  size  and  different  shapes 

•  Window  function  fixes  time  and  frequency  resolution 
-  Uncertainty  principle  dictates  ( A<)(  A/)  >  ^ 

•  Useful  for  analysis  of  narrowband  processes 

_ MITRE _ _ 


Gabor  Wavelets 

To  improve  time-localization  over  the  Fourier  transform,  one  can  multiply  the  signal  by  a 
“window”  function  first.  Since  the  Gaussian  function  has  optimal  time-frequency  spread 
according  to  the  Uncertainty  Principle  (equality  is  achieved),  Gabor  chose  the  Gaussian  as  his 
window  function.  Thus,  the  Gaussian  function,  translated  and  modulated  integral  amounts,  is 
the  complete  system  of  functions  called  the  Gabor  wavelets.  Note  that  unlike  the  Fourier 
basis,  this  system  is  not  orthogonal. 
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[  Gabor  Wavelets  (Continued)  1 


•  Basic  window  function  g,  e.g. 

g(t)  =  X[_i/2,i/2)(0  Haar 

$(<)  =  it-1/4  e-t  /2  Gabor 

•  Energy  distribution  in  time:  (Af.)2  =  II*p(0II2/IIpII2 

Energy  distribution  in  frequency:  (A/)2  =  ||/<K/)ll2/ll9ll2 

•  Window  functions  are  defined  by 

9a, b^)  —  e2*'atg(t-b) 

»„,*(/)  =  «-**/-•>*$(/  -  a) 

•  If  energy  in  time-frequency  plane  for  g  is  centered  at 
(ao,bo)  then  energy  in  time-frequency  plane  for  ga  b  is 
centered  at  (a0  4-  a,  b0  4-  b) 

i  MfTKE 


Gabor  Wavelets  (Continued) 

The  uncertainty,  or  spread,  of  a  function  in  time  and  frequency  is  defined  as  above. 
Regardless  of  which  function  is  chosen,  the  product  of  the  two  spreads  is  bounded  below 
(Heisenberg).  Furthermore,  translation  and  modulation  do  not  change  the  characteristics  of 
spread,  but  only  shift  the  center  of  concentration  in  the  time-frequency  plane  (defined  by 

means  of  the  quantities  |/p,J/|  ). 
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Gabor  Wavelets  (Continued) 


•  Windowed  transform  Wg  maps  L2(R)  — ►  L2(R 2) 

Wgs(a,b)  =  J  s(t)gab(t)dt 

•  Reconstruction  formula 

s(<)  ~  J  f  Wgs(a,b)ga>b(t)  dadb 

•  Redundancy  in  continuous-parameter  representation 

-  Some  redundancy  can  be  good 

-  Sampled  transform  can  provide  reconstruction  of  s 

s(t )  =  V^J7s(a0m,feo«)5a0m,6on(0 

m,n 

Gabor  Wavelets  (Continued) 

For  a  choice  of  a  window  function  g  with  llgll  =  1,  the  Gabor  transform  is  given  by  an 
inner  product  of  the  signal  with  a  translated,  modulated  g.  The  Gabor  transform,  as  a 
function  of  the  translation  and  modulation  parameters,  allows  exact  reconstruction  of  the 
original  signal.  In  fact,  discrete  samples  of  the  transform  can  also  provide  reconstruction,  for 
special  choice  of  lattice  parameters  oq,  bo. 
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Gabor  Wavelets  (Continued) 


\ 


•  Reconstruction  depends  on  choice  of  parameters  ao,  bo 

9m,n(t)  =  e2nta°mtg(t  -  b0n) 


MM>  I 


-j.  ■  -  ■  -  - 

•  2  4  *  t  M 


(*0,b0)-pUno 


a0bo  <  1  Frames 
a0bo  =  1  Bases  Possible 
a0b0  >  1  Undersampling 


v 
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Gabor  Wavelets  (Continued) 

For  flobo  ^  1*  reconstruction  is  possible  (g  has  to  be  well-chosen  also).  The  case  ao&o  =  1 
is  special  in  that  actual  bases  are  possible,  not  just  frames  (which  often  have  more  vectors 
than  necessary).  For  oofeo  >  1»  reconstruction  is  impossible  for  any  g. 
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Gabor  Wavelets  (Concluded) 


•  Uncertainty  Principle  dicates  that  (At)(Af)  >  — 

{0m, n}  basis  (exact  frame)  =>  ||tg(t )llll/s(/)H  =  o° 

•  Thus  good  time-frequency  localization  is  not  possible  for 
Gabor  wavelet  bases  (exact  frames) 

•  Equivalently,  good  time-frequency  localization  requires 
redundancy  in  the  Gabor  transform 


MT7RE 


Gabor  Wavelets  (Concluded) 

Unfortunately,  for  such  a  Gabor  frame  to  be  a  basis,  the  window  function  must  have  poor 
time-frequency  localization  properties — either  the  time  or  frequency  spread  diverges.  In 
particular,  the  Gaussian  window  does  not  give  rise  to  any  basis,  since  its  time  and  frequency 
spreads  are  finite  (in  fact  optimal).  In  general,  good  time-frequency  localization  with  Gabor 
wavelets  can  be  achieved  only  by  using  frames. 
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Affine  Wavelets 


•  Window  functions  are  obtained  by  dilating  and 
translating  a  single  function 

-  Functions  have  different  sizes  and  same  shape 

•  Provides  for  multiscale  analysis 

-  Good  time  resolution  at  high  frequencies 

—  Good  frequency  resolution  at  low  frequencies 

MITRE 


Affine  Wavelets 

Affine  wavelets  are  constructed  by  translations  and  dilations  of  a  single  functions.  These 
give  rise  to  a  different  tiling  of  the  time-frequency  plane  when  considering  discrete  bases.  In 
essence,  affine  wavelets  achieve  good  time  resolution  at  high  frequencies,  and  good 
frequency  resolution  at  low  frequencies  (at  the  expense  of  good  frequency  or  time  resolution, 
respectively).  For  numerous  applications,  this  choice  in  the  time-frequency  trade-off  is 
advantageous. 
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Affine  Wavelets  (Continued) 


•  General  Theory 

€  L2(R)  fixed;  set  —  a}l2rl>(at  —  b) 

:  L2{R)  -  L2(R2) 

U^s(a,6)  =  f  dt 


•  For  the  reconstruction  of  s 

s(t)  =  C^1  JJ  Wj,s(a,b)1>a,b(t)^ 

we  need  an  admissibility  condition 

c^  =  jif\<Kf)\2/\f\<V<<x> 

so  that  V>( 0)  =  0 
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Affine  Wavelets  (Continued) 

The  affine  wavelet  transform  is  given  by  an  inner  product  of  the  signal  with  a  translated, 
dilated  wavelet  It  also  allows  reconstruction  of  the  original  signal,  under  mild  conditions  on 
the  wavelet  function  (it  has  zero  mean). 


I  Affine  Wavelets  (Continued) 


•  Discrete  Version:  i/>m,n(t)  =  -  60n) 

•  Now 

so  that  the  parameter  m  determines  the  resolution  in  time 
and  frequency  while  the  parameter  n  locates  the  center 
of  energy  within  each  scale 

m  »  0  :  good  time  resolution 
m  <  0  :  good  frequency  resolution 
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Affine  Wavelets  (Continued) 

Again,  discrete  samples  of  the  wavelet  transform  can  suffice  for  reconstruction.  The 
dilation  parameter  effectively  controls  the  time-frequency  characteristics  of  the  wavelet. 


30 


Affine  Wavelets  (Concluded) 


•  Consider  the  Haar  example:  a0  =  2,  60  =  1 

4’  =  4’ 0,0  =  *[0,l/2)  ~  *[1/2,1) 


V 
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Affine  Wavelets  (Concluded) 

As  a  basic  example,  the  Haar  wavelet  and  certain  of  its  translated  and  dilated  cousins  are 
depicted. 


Wavelet  Transform 


Wj,8(a,b)  =  /  s(t)rl>a  b(t)  dt 

J—oo 


Gabor 

a  =  m,  b  —  n 
V'm.nCO  =  c2ff*mt  g(t  -  n) 


Affine 

a  =  2m,  b  =  n 
xPm,n(t)  =  2m/2rP(2mt-n) 


«(<)  =  51  W^s(m,n)^m,n(0 


n,m=— oo 


wavelet  transform 


Wavelet  Transform 

In  summary,  both  the  Gabor  and  Affine  wavelet  transforms  are  given  by  inner  product  of 
the  signal  against  a  two-parameter  family  of  functions.  Both  allow  discrete  sampling  to 
provide  reconstruction. 
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f  Multiresolution  Analysis  1 


A  multiresolution  analysis  for  finite  energy  signals  L2(R)  is  a 
collection  of  closed  subspaces  { Vm}  satisfying  the  properties 

•  Vm  C 

•  nVrn  —  {0} 

•  UVm  =  L2{R) 

•  s(t)  €  Vm  &  s(2t)  e  Vm+1 

•  There  exists  <p  so  that  {<p(t  -  n)}n€z  is  an  orthonormal 
basis  for  Vq 


Multiresolution  Analysis 

Essentially,  every  orthonormal  wavelet  basis  can  be  derived  from  a  construction  due  to 
Mallat  and  Meyer  called  multiresolution  analysis.  Assume  that  the  function  space  L2(R)  has 
a  sequence  of  nested  subspaces,  whose  elements  are  related  by  dilation  by  powers  of  2. 
Furthermore,  assume  that  one  such  space,  say  Vo,  has  a  function  (p  whose  integer  translates 
are  in  Vq.  and  form  an  orthonormal  basis. 
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Multiresolution  Analysis  (Continued) 


•  At  each  resolution  level  Vn+ 1  =  Vm(BWm 

•  Vrn  =  {s(2mt)  :  s  €  Vq} 

Wm  =  {s(2mt)  :seW0} 

•  We  thus  obtain  a  decomposition  of  L2(R): 

L2(R)  =  Vm  e  wm+1  ©  wm+2  ©  •  •  • 

=  ®Wm 

Here,  Vm  is  the  resolution  at  level  m  and  the  remain¬ 
der  WmJrl  ®  •  •  •  contains  the  high  frequency  components 
needed  to  generate  all  of  L2(R) 
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Multiresolution  Analysis  (Continued) 

Since  Vm  c  Vm+i,  one  can  introduce  the  orthogonal  complement  space  Wm,  so  that 
Vm+l  =  Vm©  Wm.  Since  fl  Vm  =  o,  and  UVm  =  L2(R)  one  has  that  ®  Wm  =  L2(R).  That 

mm  m 

is,  any  square  integrable  function  can  be  written  as  a  sum  of  components  that  live  in  the  Wm. 
In  fact,  the  spaces  Wm  are  spanned  by  the  translates  of  a  fixed  dilate  of  a  single  function,  yr. 
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Multiresolution  Analysis  (Concluded) 


•  Since  Vq  c  Vi  there  are  coefficients  {cn}  so  that 

*p(t)  =  J2Cn^(2t  -  n) 

n 

•  A  basis  for  W0  is  given  by  {ip(t  -  n)}neZ  where  ip  is  given 
explicitly  in  terms  of  <p 

VKO  =  5Z(-l)nc1_„v>(2<  -  n) 

n 

•  <p,  ip  satisfy  fv=  1  and  /  ip  =  o 


V  MfTKE 


Multiresolution  Analysis  (Concluded) 

Given  this  setup,  the  function  <p  satisfies  a  recursion  relation  called  a  dilation  equation. 
The  coefficients,  or  taps,  in  the  relation  play  an  important  role.  A  formula  exists  for 
constructing  a  function  y  from  <p  and  its  cousins,  such  that  translates  and  dilates  of  yf  form 
an  orthonormal  basis,  y  automatically  satisfies  the  requirement  of  zero  mean. 
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Example:  Daubechies  4-Tap  Wavelet 

t 

ftS 

• 

/V 

1 

*5  1  M  2  IS 

Scaling  Function 

i  e«  i  ij  >  is 

Wavelet 

l  +  y/3  3+V3 

3-V3  l-y/3 

C0  4  Cl“  4 

C2  4  C3  4 

MTTRE 

Example:  Daubechies  4-Tap  Wavelet 

The  Haar  function  satisfies  a  2-tap  dilation  equation.  Daubechies  (75  years  later) 
discovered  a  function  that  satisfies  a  4-tap  relation.  The  functions  <p,  yr  (called  scaling 
function  and  wavelet,  respectively)  in  this  case  are  very  unusual  in  that  they  are  continuous 
but  not  differentiable. 
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Example:  Daubechies  10-Tap  Wavelet 

Allowing  more  taps  in  the  dilation  equation  permits  smoother  solutions,  though  still  not 
infinitely  differentiable. 
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Example:  Non-Continuous  4-Tap  Wavelet 


Example:  Non-Continuous  4-Tap  Wavelet 

At  the  other  extreme,  even  non-continuous  wavelets  exist  (the  Haar  function  is  a  simple 
2-tap  example). 
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Fast  Wavelet  Transform 

*0 

G-@— 

yi  _ 

- 

- G  — (i2) - 

V2 

—  H-t&r2- 

— [U-0— 

V3 

- tf-0-53. 

•  ’  * 

N 

N 

N 

N 

2 

4 

8 

•  G  highpass  filter  {(-l)”^. 

-n/V2) 

H  lowpass  filter  {cn/\/ 2} 

•  Perfect  reconstruction  of  x0  from  xp,y1,y2,...,yp  with  1 

order  N  algorithm 

MfTRE 

_ ) 

Fast  Wavelet  Transform 

The  coefficients  in  the  relations  for  <p,  if/  give  rise  to  FIR  filters  H,  G,  which  serve  as 
lowpass  and  highpass  filters,  respectively.  A  sampled  signal  of  length  N  is  high  and  lowpass 
filtered  and  subsampled.  The  resulting  lowpass  signal,  of  length  N/ 2,  is  again  high  and 
lowpass  filtered,  and  so  on,  until  the  process  stops.  The  resulting  set  of  highpass  signals,  plus 
one  lowpass  term,  allow  exact  reconstruction  of  the  original  signal.  The  whole  process  is 

om 
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Fast  Wavelet  Transform  (Concluded) 


•  Analysis 

Data 

LP  Filter  Coefficients 
HP  Filter  Coefficients 
Smoothed  Data 
Detail  Data 


*oM 

h[n]  =  c[n]/y/ 2 
<?[n]  =  (-l)Ml  -  n]/V2 
=  EjfcM2*  -  *]*olfe] 
VI M  =  LfcSl 2n  -  k)x0[k] 


•  Synthesis 


where 


*oM  =  ]C*ilk]j»[fc]  +  yi[fc]y[fc] 

k 


&(/:]  =  h[2k  —  n]  and  g[k]  =  g[2k  —  n] 


/ 


Fast  Wavelet  Transform  (Concluded) 

This  slide  gives  a  synopsis  of  the  digital  wavelet  analysis  and  reconstruction  formulas. 
Note  that  the  highpass  filter  is  directly  constructed  from  die  lowpass  one. 
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Daubechies’  List 


In  1988,  Daubechies  published  a  list  of  useful  wavelet  lowpass  filters,  which  she  labeled 
as  N^n  for  integers  N  (unrelated  to  signal  length;  2^  is  the  number  of  taps).  Here  we  present 
the  4, 6, 8, 10, 12,  and  14  tap  filters. 


n 

Nhn 

N  =  2 

0 

.4829629131445341 

1 

.8365163037378077 

2 

.2241438680420134 

3 

-.1294095225512603 

N  =  3 

0 

.3326705529500825 

1 

.8068915093110924 

2 

.4598775021184914 

3 

-.1350110200102546 

4 

-.0854412738820267 

5 

.0352262918857095 

N  =  4 

0 

.2303778133088964 

1 

.7148465705529154 

2 

.6308807679398587 

3 

-.0279837694168599 

4 

-.1870348117190931 

5 

.0308413818355607 

6 

.0328830116668852 

7 

-.0105974017850690 

N  =  5 

0 

.1601023979741929 

1 

.6038292697971895 

2 

.7243085284377726 

3 

.1384281459013203 

4 

-.2422948870663823 

5 

-.0322448695846381 

6 

.0775714938400459 

7 

-.0062414902127983 

8 

-.0125807519990820 

9 

.0033357252854738 
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Daubechies’  list  2 


The  graphs  of  some  of  the  scaling  and  wavelet  functions  from  this  list  are  depicted  here. 
Note  that  the  functions  are  increasingly  smoother  down  the  list. 
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f  Time  Frequency  Localization  via  the  | 
I  Weyl  Correspondence  1 


*  Technique  for  focalizing  a  signal  in  a  region  of  the  time- frequency  (phase)  plane 

*  Generaltzee  bandpass  and  time  pass  filters  to  time-varying  bandpass  filters 

*  Depends  on  a  time- frequency  Joint  distribution  and  a  region  in  the  phase  plane 


Time  Frequency  Localization  via  the  Weyl  Correspondence 

We  now  consider  in  more  detail  the  notion  of  time-frequency  localization  of  signals.  The 
Fourier  and  wavelet  bases  gave  rise  to  time-frequency  tilings  with  a  certain  structure. 
Projecting  a  signal  onto  the  span  of  a  subset  of  these  basis  vectors  would  amount  to 
projecting  onto  a  subset  of  the  time-frequency  plane.  While  the  subsets  involved  had  a 
rectangular  structure  previously,  we  will  now  study  such  projection  operators  that  relate  to 
arbitrary  domains. 

To  begin  with,  projecting  onto  a  subset  of  the  time-frequency  plane  is  a  common 
operation  in  filtering  theory.  For  example,  a  bandpass  filter  projects  onto  a  strip  domain 
parallel  to  the  time  axis.  Similarly,  timepass  filtering  (restricting  the  signal  to  a  time  domain) 
can  be  viewed  as  projecting  onto  a  strip  parallel  to  the  frequency  axis.  But  how  can  one 
project  onto  an  arbitrary  domain?  And  which  functions  are  well-localized  in  such  domains? 
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Preliminaries 


Sampling  Theorem 

f  (x)  s  Lf  (n)  sine  (  x  (x  -  n)) 
lor  f  (x)  bandllmlted  to  [-1/2, 1/2] 


w 

10 

i 

tpfcfis 

. 

V.0.V.,.V.V.'.X.0.V.V.V.V 

•10 

Heisenberg  Uncertainty 

II  X f  ||  ||  Df||  £1/4*  ||f  ||2 
D=  1/2 ici  d/dx 


Extremals: 

1/4  2 

$(X)s2  «Cp(-  KX  ) 

exp(2jdax)  $(r(x-b)) 


V 


MURE 


Preliminaries 

As  a  preliminary,  consider  bandlimited  function.  The  sine  function  is  bandlimited.  In 
fact,  its  integer  translates  form  an  orthonormal  basis  for  bandlimited  functions.  Even  more, 
the  coefficients  involved  in  expanding  any  bandlimited  function  in  this  basis  are  actually  the 
values  of  the  function  at  the  integer  points  (Sampling  Theorem).  Note  that  one  has  one  basis 
function  per  unit  area  in  the  band  domain-this  is  the  Nyquist  density  of  states. 

The  bandlimited  sine  functions  are  analytic  functions,  making  them  very  special. 
However,  though  limited  in  frequency,  they  have  very  slow  decay  in  time  (1/x).  Functions 
that  are  well-localized  in  both  time  and  frequency  (e.g.,  the  product  of  their  time  and 
frequency  spreads  is  small)  tend  to  be  even  better.  For  example,  the  extremals  of  the 
Uncertainty  Principle  are  the  (modulated,  dilated,  and  translated)  Gaussians,  which  have 
exponential  decay  in  time  and  frequency,  in  addition  to  being  analytic. 
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Wavelet  Transforms 

•  Affine  wavelet  transform 

Tf(a,b)  =  J  f(x)  fa  g(a  x  -  b)  dx 

•  Gabor  (W-H  wavelet)  transform 

• 

• 

• 

• 

• 

• 

t 

• 

b 

Gf(a,b)  =  J  f(x)  exp(2itiax)  g(x  -  b)  dx 

• 

t 

MITRE 

• 

§  b 

• 

Wavelet  Transforms 

For  square-integrable  signals  that  are  not  necessarily  bandlimited,  there  is  still  a  discrete 
representation  available,  now  in  the  Affine  or  Gabor  Wavelet  Transform.  In  this  case,  the 
sampling  is  on  a  two-dimensional  lattice  (the  affine  case  is  also  on  a  lattice,  when  log  a  is 
used  as  a  parameter  instead  of  a). 
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Time-Frequency  Decomposition 

In  terms  of  time-frequency  concentration,  the  Fourier  and  wavelet  bases  have  rectangular 
domains  of  support.  In  the  Fourier  case,  the  rectangles  are  of  identical  shape,  while  in  the 
wavelet  case  they  vary,  though  still  of  the  same  area.  To  find  basis  functions  that  have  more 
general  support  domains,  we  study  a  class  of  time-frequency  localization  operators. 
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Prolate  Spheroidal  Wavefunctions 


Landau,  Pollack  and  Sleplan  (1961) 
For  f  with  ||f||  =  1, 

max  ||ABf|| 

Equivalently,  maximize: 

(ABf,  ABf)  *  (ABf,  Bf  > 

=  <BABf,  f  > 


•  Thus,  seek  maximum  eigenvalue  (and  eigenfunction)  of 
the  positive,  self-adjoint  operator  BAB. 

•  Succeeding  eigenvalues  are  maximal  in  the  space 
orthogonal  to  the  previous  eigenfunctions.  Get  the 
well  known  prolate  spheroidal  wavefunctions. 
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Prolate  Spheroidal  Wavefunctions 

In  1961,  Landau,  Pollack,  and  Slepian,  at  Bell  Labs.,  defined  a  localization  operator  for  a 
rectangular  region  in  time-frequency  space.  Their  operator  was  just  a  concatenation  of  the 
bandpass,  timepa'",  and  bandpass  operators.  It  emerged  that  their  operator  commuted  with  a 
Laplace  type  operator,  whose  solutions  were  well-known:  the  prolate  spheroidal 
wavefunctions. 
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Wigner  Distribution 


Tlme-freq.  shift:  p(a,b)f(x)  =  exp(7riab)exp(27dax)  f(x+b) 
Ambiguity  Fn.:  A(f,g)(a,b)  =  <  p(a,b)f,g> 

Wigner  Dist.:  W(f,g)  =  A(f,g)* 

energy  representation:  satisfaction  of 


l«q 

lW 

1 

tin 

|t(x)2 

1 
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Wigner  Distribution 

To  build  an  operator  corresponding  to  more  general  domains,  we  have  to  proceed  in  a 
different  way.  First,  consider  the  Wigner  distribution,  which  is  nothing  else  than  the  Fourier 
transform  of  the  Ambiguity  function.  The  Wigner  distribution  Wigner  distribution  of  a 
signal  is  essentially  a  distribution  of  the  signal’s  energy  in  the  time-frequency  plane.  The 
integral  of  the  Wigner  distribution  over  the  plane  yields  the  energy  of  the  signal,  while 
integrals  along  lines  parallel  to  the  axes  yield  the  power  and  the  power  spectral  density, 
respectively. 
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Weyl  Correspondence 


Given  a  symbol  function  S(w,t), 
define  an  operator  Lg  by 


<l-sf.9>=  J  S(w,t)  W(f,g)(w,t)dwdt 

~i  Q 

Example:  S  =  X 


f 

General 


Wcyl  Correspondence 

Using  the  Wigner  distribution,  H.  Weyl  defined  an  operator  as  above  for  any  “symbol” 
function  S.  We  are  interested  in  the  case  when  S  is  the  characteristic  function  of  a  set  Q  (i.e., 
S(w,  t)  *  1  if  (w,  t)  e  Q,  0  else). 
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Cohen's  Class  Operators 

Cohen's  class  distribution: 

P(f,g)(w,t)=  J  W(f,g)(a,b)  s(w -a,t-b)dadb 

p 

Cohen's  class  operator  Lg: 

<L^f,  g>=  J  S (w , t)P (f , g )( w , t )dwdt 

*-S  =  S'  S(w,t)=  J  S(w,t)  s.(a-w,b-t)dadb 

*_(a,b)  *  »{-a,  -b) 

,  MTTFE 


Cohen’s  Qass  Operators 

In  fact,  such  a  procedure  for  defining  operators  can  be  generalized  from  the  Wigner 
distribution  to  any  one  of  the  so-called  Cohen’s  class  of  time-frequency  distributions. 
However,  such  a  generalization  is  illusory,  and  merely  amounts  to  a  change  of  the  symbol 
function.  The  utility  of  working  with  distributions  other  than  the  Wigner  distribution  is 
mainly  in  the  interpretation. 
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Daubechies  Operator 

In  1988,  Daubechies  defined  such  a  localization  operator,  in  terms  of  the  “spectrogram” 
distribution,  which  is  the  modulus-squared  of  the  short-time  Fourier  transform  with  a 
Gaussian  window.  The  spectrogram  is  a  well-known  tool  in  speech  processing. 
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Weyl  vs.  Daubechies 

In  the  case  of  localization  onto  a  circular  disc,  it  turns  out  that  the  eigenfunctions  for  both 
the  Weyl  and  Daubechies  operators  are  identical:  the  Hermite  functions.  In  this  case,  one 
can  compare  the  eigenvalues.  While  there  is  a  marked  drop-off  near  100  (the  area)  in  both 
cases,  there  arc  two  differences.  The  drop-off  is  faster  in  the  Weyl  case,  but  there  is  a  Gibbs- 
like  ringing  near  the  zero  eigenvalue.  Note  that  for  numerical  purposes,  the  Weyl  operator 

allows  a  sharper  focus  of  energy,  if  projecting  onto  eigenvectors  with  eigenvalues  >  -. 

2 
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f  Prolate  Spheroidal  Wavefunctions  1 


Localization  on  a  4x2 
rectangle 

Eigenfunctions  = 
Prolate  spheroidal 
wavefunctions 
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Prolate  Spheroidal  Wavefunctions 

The  next  three  slides  consider  projecting  onto  a  rectangle,  of  area  =  8.  Here  the  first  nine 
prolate  spheroidal  wavefunctions  are  presented.  Because  these  functions  are  exactly 
bandlimited,  they  have  a  slow  drop-off  in  time  (e.g.,  \/x),  as  evidenced  by  the  last  three 
functions. 
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Weyl  Eigenfunctions 


\ 


Localization  on  a  4x2 
rectangle 

Eigenfunctions  = 
New  Weyl  solutions 


First  9  Weyl  eigenfunctions 
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Weyl  Eigenfunctions 

When  using  the  Weyl  correspondence,  one  finds  eigenfunctions  with  a  very  similar 
appearance  to  the  prolate  functions,  although  the  drop-off  in  time  is  better.  MITRE  research 
mathematicians  were  able  to  derive  bounds  on  the  rate  of  decay  of  the  solutions. 
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Daubechies  Eigenfunctions 


Localization  on  a  4x2 
rectangle 

Eigenfunctions  = 
New  Daubechies 
solutions 


First  9  eigenfunctions 


Eigenvalues 


I  if  IS 


MITRE 


Daubechies  Eigenfunctions 

In  the  Daubechies  case,  the  solutions  appear  similar  (again),  but  with  even  better  decay. 
Note  that  by  contrast,  the  eigenvalues  decay  slowly,  indicating  that  the  total  energy  is  spread 
to  many  modes. 


58 


New  Results  (Wigner  Distribution  Case)  1 


•  For  localization  onto  a  general  bounded  domain  with  a 
piecewise  (once)  differentiable  boundary,  we  prove 
that: 


-  Wigner  solutions  are  analytic 

-  They  have  exponential  decay,  for  any  positive  exponent 

-  Their  Fourier  transforms  satisfy  identical  conditions 

•3/4 

-  Eigenvalues  decay  at  the  rate  of  n 
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New  Results  (Wigner  Distribution  Case) 

The  qualitative  observations  made  earlier  are  confirmed  by  theorems,  proved  in 
collaboration  with  Jayakumar  Ramanathan  of  E.  Michigan  University,  that  the  eigenfunctions 
of  the  Weyl  operator  are  analytic  and  have  exponential  decay,  under  mild  hypotheses 
regarding  the  domain  of  localization.  [J.  Ramanathan  and  P.  Topiwala,  ‘Time-Frequency 
Localization  and  the  Weyl  Correspondence,”  SIAM  Journal  of  Mathematical  Analysis, 
September  1993.] 
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New  Results  (Spectrogram  Case) 

•  For  localization  onto  a  general  bounded  measurable 
domain,  we  prove  that: 

-  Spectrogram  solutions  are  analytic  and  of  quadratic  exp.  decay 
•  They  and  their  derivatives  decay  faster  than  l/|x|  to  any  power 

-  Their  Fourier  transforms  satisfy  identical  conditions 

-  Eigenvalues  decay  faster  than  1/n  to  any  power 

"  The  number  of  eigenvalues  above  a  fixed  number  y,  0<y<l, 

Is  asymptotically  equal  to  the  area  of  the  domain,  upon 
rescaling  the  domain  -  Nyquist  density. 

MITRE 


New  Results  (Spectrogram  Case) 

For  the  spectrogram  case  (Daubechies  operator),  the  results  are  even  stronger:  the 
solutions  are  analytic  and  have  quadratic  exponential  decay,  under  very  general  conditions  on 
the  domain.  We  also  have  a  result  along  the  lines  of  Nyquist  density:  for  a  given  domain, 
the  number  of  eigenvalues  greater  than  a  fixed  number  y,  0  <  y  <  1,  is  asymptotically  equal  to 
the  area,  as  the  domain  is  dilated  radially.  These  sets  of  results  confirm  the  intuition  that 
signals  that  are  well-localized  in  the  time-frequency  plane  satisfy  strong  regularity  and  decay 
conditions.  Furthermore,  once  found,  they  lead  to  simple  projection  operators,  by  projecting 

onto  the  span  of  the  eigenfunctions  with  eigenvalues  greater  than  some  cutoff  (e.g.,  -) 

[J.  Ramanathan  and  P.  Topiwala,  “Time-Frequency  Localization  and  the  Gabor 
Transformation,”  Journal  of  Applied  and  Computational  Harmonic  Analysis  (to  appear).] 
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Applications 


•  1-D  Signal  Filtering 

•  1-D  Transient  Detection 

•  Image  Compression 
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Applications 

Wc  now  consider  some  applications  of  the  ideas  and  tools  developed  so  far.  We  will 
consider  only  three  topics:  noise  and  interference  filtering,  transient  detection,  and  image 
compression. 
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Noise  Filtering 


Example  of  noise  filtering  from  a 
linear  FM  (chirp)  signal 


rm*mt 


« 


Wigner  distribution  of  a)  chirp 
signal;  b)  chirp  in  0  dB  noise; 
c)  Weyl  filtered  noisy  signal 
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Noise  Filtering 

In  this  example,  the  Wigner  distribution  of  a  chip  signal  is  given  in  (a),  (b)  is  the  Wigner 
distribution  of  the  chirp  with  0  dB  Gaussian  noise  added.  Note  that  while  the  Wigner 
distribution  of  the  chirp  is  a  straight  line,  corresponding  to  the  slope  of  the  chirp,  the  noise 
has  energy  distributed  all  over  the  time-frequency  plane,  (c)  is  the  Wigner  distribution  of  the 
noisy  signal  after  it  has  been  Weyl  filtered,  using  a  narrow  strip  along  the  chirp. 
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Noise  Filtering  (cont'd) 


Plots  of  signals 


*  w 


a)  chirp  signal; 

b)  chirp  In  0  dB  noise; 

c)  Weyl  filtered  noisy 
signal 


Noise  Filtering  (Continued) 

This  slide  shows  the  actual  signals.  Note  that  while  (c)  is  not  an  exact  reconstruction  of 
(a),  it  is  a  remarkable  improvement  on  (b). 


L 


63 


Noise  Filtering  (final) 


chirp  signal 


noisy  chirp 


filtered  noisy  signal 


•  Enhanced  Wigner  distribution  Intensity  plots 
using  windowing  of  signals  (cos2). 
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Noise  Filtering  (Final) 

The  ringing  effect  in  the  Wigner  plots  can  be  eliminated  by  simple  windowing  of  the 
signal-here  a  cos2*  window  is  used. 
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1-D  FT  Signal,  R=1 

In  this  example,  we  will  consider  detection  of  a  high-frequency  low  magnitude  transient 
signal  in  a  background  smooth  process.  Here  a  Gaussian  pulse  is  used  as  background,  and  a 
high-frequency  sinusoid  as  the  interesting  transient.  We  will  compare  the  effects  of  using 
FFT  vs.  wavelet  processing,  in  particular  compression,  on  the  detectability  of  the  transient 
In  this  slide,  the  original  signal  is  reconstructed  exactly  under  FFT  processing  (FFT  and  IFFT 
applied),  with  compression  factor  #=1  (no  compression).  Note  that  compression  is 
considered  here  as  a  method  of  feature  extraction. 
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1-D  FT  Decomposition,  R=1 


1-D  FT  Decomposition,  R=1 

Here  the  real  and  imaginary  parts  of  the  Fourier  transform  (FFT)  are  displayed.  The 
high  frequency  transient  is  coded  mainly  near  the  middle  of  these  plots. 
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1-D  WT  Signal,  R=1 


Similarly,  the  Wavelet  transform  (WT)  reconstructs  the  signal  exactly  without 
compression. 
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1-D  WT  Decomposition,  R=1 


Here,  the  first  two  lowpass  and  highpass  components  are  depicted  (out  of  5  total 
computed).  Note  that  unlike  the  FFT,  the  character  of  the  transient  and  Gaussian  pulse  are 
partly  reflected  in  these  frames.  In  particular,  the  time-domain  on  which  the  transient  occurs 
is  clearly  visible  in  the  first  highpass  component  In  the  following  WT  slides,  attention 
should  be  focused  on  that  region.  Note  that  the  lowpass  filter  acts  as  a  local  integrator,  which 
is  why  the  transient  is  not  detected  in  that  component;  the  highpass  filter  is  a  differentiator, 
which  reacts  strongly  to  the  transient. 
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1-D  FT  Signal,  R=10 


Here  10:1  compression  has  been  applied  in  the  Fourier  domain  before  inversion.  All 
traces  of  the  transient  are  lost;  the  resulting  signal  is  smooth  but  significantly  deformed  (in  a 
mean-square  error  sense). 
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1-D  FT  Decomposition,  R=10 

1 

Had  pvt  of  Fourier  transform 
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1-D  FT  Decomposition,  R=10 

Compression  was  implemented  by  setting  the  lowest  90%  of  the  terms  (in  magnitude)  to 
0.  This  had  the  effect  of  deleting  all  the  high-frequency  coefficients  (which  would  normally 
reside  in  the  center  of  these  plots  due  to  the  cyclic  nature  of  the  FFT). 
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1-D  WT  Signal,  R=10 

By  contrast,  the  wavelet  transform  retains  some  features  of  the  transient  under  10:1 
compression.  In  general,  the  signal  is  closer  to  the  original  in  a  mean-square  sense  than  the 
Fourier  reconstructed  signal. 
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1-D  WT  Decomposition,  R=10 

In  the  wavelet  decomposition,  the  10:1  compression  still  retains  the  effect  of  the  transient 
in  the  first  highpass  component.  Thus,  a  successful  detector  for  this  can  be  designed  based 
on  the  first  highpass  component  by  setting  a  threshold.  Again,  the  time  of  occurrence  of  the 
transient  is  directly  available,  as  well  as  a  certain  measure  of  its  intensity.  We  note  that  other 
types  of  transients  may  be  effectively  represented  in  other  components  of  the  wavelet 
decomposition,  depending  on  their  nature. 
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1-D  FT  Signal,  -3dB  N,  R=10 


1-D  FT  Signal,  -3dB  N,  R=10 

We  continue  the  experiment  by  adding  -3dB  Gaussian  noise,  which  has  no  effect  in  the 
F'FT  case. 
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1-D  FT  Decomposition,  -3dB  N,  R=10 


1-D  FT  Decomposition,  -3dB  N,  R=10 

The  noise  is  filtered  away  along  with  the  transient  under  compression  in  the  Fourier 
domain. 
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1-D  WT  Signal,  -3dB  N,  R=10 

By  contrast,  the  WT  continues  to  retain  features  of  the  transient,  while  filtering  most  of 
the  noise. 
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1-D  WT  Decomposition,  -3dB  N,  R=10 

Again,  a  threshold-based  transient  detector,  using  the  first  highpass  component,  would  be 
successful  here. 
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Lena  Original 


We  now  consider  image  compression,  which  is  perhaps  the  most  promising  application  of 
wavelets  to  date.  Here  an  original  512x512  pixel  8-bit  grayscale  image  (“Lena”)  is 
presented.  This  image  is  an  industry  standard  in  image  compression.  It  contains  a  variety  of 
textures,  from  the  smooth  skin  to  the  ruffled  feathers  in  the  hat.  “Lena”  and  “House”  are 
public  domain  images,  obtained  from  a  database  at  the  University  of  Sourthem  California. 
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LENA  ORIGINAL 
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Lena  JPEG  44:1 


This  slide  presents  a  highly  compressed  version  of  the  Lena  image,  using  an  algorithm 
that  is  now  the  international  standard  for  still  image  compression-the  Joint  Photographic 
Experts  Group  (JPEG)  algorithm.  We  used  public  domain  JPEG  software  courtesy  of  the 
Independent  JPEG  Group.  This  algorithm  divides  the  image  into  8x8  pixel  blocks,  computes 
a  Discrete  Cosine  Transform  (DCT)  on  each  block,  and  then  requantizes  as  necessary  to 
achieve  the  desired  compression  ratio.  Huffman  coding  is  also  used  to  achieve  the  final 
entropy  coding  of  the  quantized  blocks.  At  44:1  compression,  the  image  displays 
considerable  tiling  as  many  8x8  blocks  are  reduced  to  a  single  grayscale.  Note  that  the 
texture  on  the  face,  shoulders,  and  woodwork  in  the  background  is  lost. 

Note  that  blocking  of  the  image  data  is  necessary  for  DCT-based  algorithms  to  (a) 
preserve  locality  of  image  data  under  the  transform,  and  (b)  bound  the  computational 
complexity,  which  is  0(N  log  N}- blocking  bounds  the  log  N  factor. 
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LENA  JPEG  44:1  2/28/93 
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Lena  Wavelet  45:1 


At  45:1,  the  Wavelet  Transform  (WT)  retains  more  of  the  structure  of  the  original, 
including  edge  as  well  as  texture  information.  In  general,  the  WT  tends  to  retain  hard  edges 
best,  and  will  sacrifice  textural  information  at  high  compression.  Note  that  blocking  of  the 
image  data  is  completely  unnecessary  in  this  case,  since  the  WT  preserves  locality  and  has 
O(N)  complexity.  We  used  public  domain  wavelet  software  courtesy  of  the  MIT  Media 
Laboratory. 


81 


LENA  WAVELET  45:1  2/28/93 


Lena  JPEG  76:1  (Maximum) 


At  76:1,  the  maximum  compression  possible  on  this  image  under  JPEG,  the  image  is 
unrecognizable.  In  fact,  it  could  pose  as  excellent  cubist  art. 
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LENA  JPEG  76:1  (MAXIMUM)  2/28/93 


Lena  Wavelet  85:1 


By  contrast,  at  the  extremely  high  compression  of  85:1,  the  image  is  still  recognizable. 
The  eyes  are  still  well-preserved,  which  is  important  for  recognition.  Among  the  artifacts 
created  by  the  WT  are  aliasing  effects  around  edges. 
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LENA  WAVELET  85:1  2/28/93 
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House  Original 

As  a  second  example,  we  treat  an  aerial  view  of  a  residential  area.  (Unfortunately,  this 
reproduction  is  very  dark.) 
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HOUSE  ORIGINAL  2/28/93 


House  JPEG  40:1 


At  40:1,  JPEG  produces  a  low-quality  blocky  image.  The  sidewalks  are  still  visible,  but 
the  roads,  trees,  and  many  houses  are  obscured  and  tiled.  The  tiling  is  especially  awkward 
for  automatic  target  recognition  purposes,  since  it  introduces  artificial  straight  lines. 
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House  Wavelet  39:1 


At  39:1,  the  wavelet  image  is  also  of  poor  quality,  indicating  that  the  advantages  of 
wavelet  compression  over  JPEG  are  not  uniform  over  all  image  types.  Again,  the  sidewalks 
are  better  preserved  then  the  streets,  trees.  Straight  edges  continue  to  be  associated  to  houses, 
though. 
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HOUSE  WAVELET  39:1  2/28/93 
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Wavelet  Code 


We  conclude  with  some  simple  MATLAB  code  for  wavelet  processing  of  1-D  signals, 
including  compression.  This  code  uses  the  Haar  wavelet,  the  simplest  of  the  wavelet  filters. 
We  have  more  sophisticated  code,  including  C  codes  for  high  performance  wavelet  signal 
processing,  for  one  and  two-dimensional  signals. 
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%  Pankaj  Tcpiwala,  The  MITRE  Corp.,  Burlington  Rd.,  Bedford,  MA  01730 
%  This  program  confutes  the  wavelet  transform  with  the  Haar  wavelet 
%  up  to  level  5;  compression  ratio  R.  Load  data  file  'data'. 

%  Data  file  must  have  length  power  of  2. 

clear; 

clg; 

%  set  ccnpression  ratio... 

R=l; 

%  retrieve  (or  create)  data... 

load  data; 

x=data'; 

%  practice  data:  t=l:128;  x=  exp(-  (abs  (t-64)  .A2) /20)  +.l;x=x'; 
NHLength  (x)  ; 

L=round(M/2)  ; 
alpha=l/sqrt  (2) ; 
hl=zeros  (L,M)  ; 
gl=zeros  (L,M) ; 

%  form  the  first  filter  matrices... 
for  i=l:L, 

hl(i,2*  (i-l)+l)=nalpha; 
gl  (i,2*  (i— 1)  +l)=alpha; 
hl(i,2*i)=alpha; 
gl  (i,2*i)  =-alpha; 
end 

%  form  the  remaining  filter  matrices  to  level  5... 

h2=hl(l:L/2,l:L) ; 

g2=gl(l:L/2,l:L) ; 

h3=h2(l:L/4,l:L/2) ; 

g3=g2(l:L/4,l:L/2); 

h4=h3(l:L/8,l:L/4)  ; 

g4=g3(l:L/8,l:L/4) ; 

h5=h4(l:L/16,l:L/8) ; 

g5=g4(l:L/16,l:L/8) ; 

%  do  the  orthonormal  wavelet  decomposition  to  order  5  . . . 

cl=hl*x; 

dl=gl*x; 

c2=h2*cl; 

d2=g2*cl; 

c3=h3*c2; 

d3=g3*c2; 
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c4=h4*c3; 

d4^j4*c3/ 

c5=h5*c4; 

d5=g5*c4; 

%  plot  the  first  two  lew  and  high  pass  versions . . . 
subplot  (221)  /plot  (cl) ; title  ( 'First  low  pass ' ) / 
subplot  (222)  ,plot  (dl) ; title  ( 'First  high  pass ' ) ; 
subplot  (223)  /  plot  (c2)  /title  ( '  Second  low  pass ' ) ; 
subplot  (224)  ,plot  (d2)  /  title  ( '  Second  high  pass ' )  / 
pause 
clgz 

%  PRESS  THE  SPACE  BAR  TO  CONTINUE. . . 

%  threshold  carpression  routine... 
yss[dl/d2/d3/d4/d5/c5]  / 
z=abs  (y)  / 
z=sort (z) / 

P=round  (M*  ( 1-1/R) )  /  if  (R=l ) ,  P=1  /  end 
val=z  (P)  / 
res=zezos  (M,  1)  / 
for  j=l:M, 

if  abs(y(j)  )>=val,  res(j)=y(j)/  else  res(j)=0/  end 
end 

%  reconstruction... 
dl=res(l:L)  / 
d2=res  (L+l  :I/fL/2)  / 
d3=res  (I^-L/2+l:7*L/4)  / 
d4=res  (7*L/4+l :  15*L/8)  / 
d5=res  (15*L/8+l:31*L/16)  / 
c5=res(31*L/16+l:M)/ 


c4=h5'*c5  +  g5'*d5/ 
c3=h4'*c4  +  g4'*d4/ 
c2=h3'*c3  +  g3'*d3/ 
cl=h2'*c2  g2'*d2/ 

xnew=hl'*cl  +  gl'*dl/ 


%  plot  the  original  against  the  reconstructed  signal . . 
subplot (211) /  plot (x) /title ('Original  signal’)/ 
subplot  (212)  /  plot(xnew)/ 
title  ( 'Reconstructed  signal ' )  / 
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